Accepted for publication in ApJ 

Preprint typeset using I^T^X style cmulatcapj v. 04/20/08 



ON 
O 
O 
(N 

X> 
IX, 



I < 

o ■ 



> 
oo 
o 
o 

(N 

(N 
O 
On 
O 



GRAIN GROWTH AND DENSITY DISTRIBUTION OF THE YOUNGEST PROTOSTELLAR SYSTEMS 
Woojin Kwon 1 , Leslie W. Looney 1 , Lee G. Mundy 2 , Hsin-Fang Chiang 1 , and Athol J. Kemball 3 

Accepted for publication in ApJ 

ABSTRACT 

We present dust opacity spectral indexes (f3) of the youngest protostellar systems (so-called Class 
sources), L1448 IRS 2, L1448 IRS 3, and L1157, obtained between the A = 1.3 mm and 2.7 mm 
continua, using the Combined Array for Research in Millimeter- wave Astronomy (CARMA). The un- 
precedented compact configuration and image fidelity of CARMA allow a better detection of the dust 
continuum emission from Class sources, with a less serious missing flux problem normally associated 
with interferometry. Through visibility-modeling at both A = 1.3 mm and 2.7 mm simultaneously, as 
well as image- and visibility-comparison, we show that /3 of the three Class sources are around or 
smaller than 1, indicating that dust grains have already significantly grown at the Class stage. In 
addition, we find a radial dependence of (3, which implies faster grain growth in the denser central 
regions and/or dust segregation. Density distributions of the Class sources are also addressed by 
visibility- modeling. 

Subject headings: circumstellar matter — stars: individual (L1448 IRS 2, L1448 IRS 3, L1157) 



1. INTRODUCTION 

Although dust grains are only about one hundredth of 
the interstellar medium by mass, they play crucial roles 
for star formation, planet formation, and furthermore the 
origin of life. They are essential places to form and store 
molecules, and they are the main ingredient to form ter- 
restrial planets, as well as playing a role in the heating 
and cooling mechanisms during star and planet forma- 
tion. 

The dust opacity 4 spectral index (/?) is related to dust 
properties. It depends on dust grain sizes, compositions, 
and shapes (e.g., Pollack et al. 1994; Draine 2006). In 
particular, it is largely sensitive to grain sizes; larger 
grains give smaller j3 (e.g., Draine 2006). Many obser- 
vational studies at infrared and millimeter wavelengths 
toward T Tauri circumstellar disks have reported smaller 
values of (3 (~ 1.0) (e.g., Andrews & Williams 2007) 
compared to that of the interstellar medium (~ 1.7) 
(Finkbeiner et al. 1999; Li & Draine 2001). In the sense 
that dust grains may develop terrestrial planets, it is very 
encouraging to see signatures of larger dust grains in T 
Tauri disks, evolved young stellar objects (YSOs), com- 
pared to grains in the interstellar medium. 

However, it is not clear when the dust grain growth 
responsible for the opacity spectral index (3 ~ 1 mainly 
occurs. For example, while Andrews & Williams (2005) 
reported grain growth along the YSO evolution from 
Class I to Class II, using spectral energy distributions 
over A = 1.3 mm and submillimeter data, Natta et al. 
(2007) did not find such a tendency (a systematic varia- 
tion of (3) . To distinguish when dust grains mainly grow 

1 Department of Astronomy, University of Illinois, 1002 
West Green Street, Urbana, IL 61801; wkwon@illinois.edu, 
lwl@illinois.edu, hchiang2@illinois.edu 

2 Department of Astronomy, University of Maryland, College 
Park, MD 20742; lgm@astro.umd.edu 

3 National Center for Supercomputing Applications, Univer- 
sity of Illinois, 1205 W. Clark Street, Urbana, IL 61801; akem- 
ball@illinois.edu 

4 Dust "emissivity" has also been used in literatures from the 
viewpoint of dust thermal "emission" . 



up to the sizes for /? ~ 1, Class YSOs are the best tar- 
gets to examine. Class YSOs are at the starting point 
of low-mass star formation and they are well defined. 
They have more massive envelopes than or comparably 
massive envelopes to their central compact objects (e.g., 
Andre et al. 1993). They are also characterized with 
well-developed bipolar outflows. Earlier stages such as 
starless cores might be another good target but they are 
hardly confined. Their physical conditions including age 
have a much larger scatter than Class sources. In ad- 
dition, they are not all expected to form stars. 

In fact, no definitive answer has been given to the opac- 
ity spectral index (3 of Class sources so far. It is an- 
other reason that this study is needed beyond the grain 
growth point of view. There are some previous studies 
about the flux density spectral indexes of Class sources, 
which are related to the dust opacity spectral indexes, 
although they have not focused on dust properties (e.g., 
Hogerheijde & Sandell 2000; Shirley et al. 2000). How- 
ever, these studies used submillimeter to 1.3 mm wave- 
lengths, which is near the range of peak intensities at 
envelope temperatures (~ 30 K), so the Rayleigh- Jeans 
approximation is invalid. In that case, the estimate of 
(3 is sensitive to the envelope temperature, which causes 
relatively large uncertainties in the (3 estimate. In ad- 
dition, optical thickness can cause another uncertainty, 
since Class YSO envelopes can be optically thick at 
submillimeter wavelengths. On the other hand, Harvey 
et al. (2003) obtained (3 ~ 0.8 toward the Class YSO 
B335 using A = 1.2 mm and 3 mm interferometric data, 
while carrying out modeling to test density distribution 
models of star formation. However, they did not have 
a good data set with comparable uv coverage at both 
wavelengths to discuss the (3 in detail. In other words, 
there are no reliable (3 estimates of Class YSOs. As 
a result, many studies to estimate masses from spectral 
energy distributions (SEDs) and/or to constrain density 
distributions have assumed (3^1 (e.g., Looney et al. 
2003) or considered a possible range of (3 (e.g., /3 = 1 — 2, 
Chandler & Richer 2000). 
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Radio intcrfcromctry at millimeter wavelengths is the 
best means to investigate the j3 of Class YSOs. As 
mentioned, optical thickness and dust temperature issues 
cause large uncertainties at shorter wavelengths. On the 
other hand, contamination of non-thermal continuum in- 
creases with wavelength, so it is not negligible at longer 
centimeter wavelengths. In addition, considering enve- 
lope sizes of Class YSOs and their environments (nor- 
mally they are within extended molecular clouds), single 
dish observations are not appropriate due to their lack 
of angular resolution and the contamination of molecular 
clouds. In contrast, interferometers provide high angular 
resolution and resolve out the emission from the large- 
scale molecular cloud. However, they may also resolve 
out emission from the Class envelopes. This is caused 
by limited uv coverage, particularly due to the absence 
of short baselines and zero-spacing. For these reasons, 
interferometers with good uv coverage are required. The 
recently commissioned Combined Array for Research in 
Millimeter- wave Astronomy (CARMA) provides the best 
opportunity with its unprecedented compact configura- 
tion and image fidelity (Woody et al. 2004). 

In this paper, we present dust opacity spectral indexes 
(3 of Class sources (L1448 IRS 2, L1448 IRS 3, and 
L1157) in order to tackle when the dust grain growth 
responsible for (3 ~ 1 mainly occurs: before or after the 
Class stage. We do a parametric modeling in uv space 
to address the (3 values, as well as image and visibility 
comparisons. In addition, we examine power-law density 
indexes via modeling. First, we discuss our observations 
and data reduction, focusing on how well our CARMA 
data incorporate with this study. Afterward, we show our 
results in images, uv visibilities, and visibility modelings. 
At the end, we discuss the implications of our results. 

2. TARGET YSOS 

We have carried out observations of three Class YSO 
regions (L1448 IRS 2, L1448 IRS 3, and LI 157) using 
CARMA in the A = 1.3 mm and 2.7 mm continuum. 
These three targets are well defined as Class YSOs 
by previous studies (e.g., Shirley et al. 2000; O 'Linger 
et al. 1999). L1448 IRS 2 and IRS 3 are located in the 
dark cloud L1448 of the Perseus molecular cloud com- 
plex at a distance of 250 pc. They were first revealed by 
IRAS observations (Bachillcr & Ccrnicharo 1986). L1448 
IRS 3 is the brightest infrared source in the dark cloud 
and has three Class sources (3A, 3B, and 3C), revealed 
by radio intcrfcrometric observations (Curiel et al. 1990; 
Terebey & Padgett 1997; Looney et al. 2000). Kwon 
et al. (2006) also studied the binary system of 3A and 
3B, the two interacting bipolar outflows, and the mag- 
netic field in the region, using polarimetric observations 
of the Berkeley Illinois Maryland Association (BIMA) 
array in the A = 1.3 mm continuum and CO J = 2 — > 1 
transition line. 

On the other hand, L1448 IRS 2 at - 3' west of IRS 
3 has not been focused on very much due to its weaker 
brightness. However, O'Linger et al. (1999) identified 
it as a Class YSO, using far-infrared up to millimeter 
continuum observations. In addition, recent deep Spitzcr 
Space Telescope (SST) IRAC observations have shown a 
large bipolar outflow spanning over 5' (Tobin et al. 2007) . 
CARMA observations in CO J = 2 1 and J = 1 -> 
transitions also show a well-developed bipolar outflow 



(Kwon et al. 2009). 

L1157 is a dark cloud in Cephcus. The distance is 
not well known but it is arguably about 250 pc (Looney 
et al. 2007). Its envelope and large bipolar outflow 
have been studied by radio single dish and interfero- 
metric observations (e.g., Bachiller et al. 2001; Gueth 
et al. 2003; Beltran et al. 2004). The bipolar outflow 
is known as chemically active, since various molecules 
have been detected and interestingly there is an abun- 
dance gradient that cannot be explained purely by ex- 
citation temperature differences (Bachiller et al. 2001). 
Recently, a flattened envelope has been detected in ab- 
sorption against polycyclic aromatic hydrocarbon (PAH) 
background emission by deep SST IRAC observations 
(Looney et al. 2007). 

3. OBSERVATIONS AND DATA REDUCTION 

We have carried out A = 1.3 mm and 2.7 mm contin- 
uum observations towards three Class sources, L1448 
IRS 2, L1448 IRS 3, and LI 157, using CARMA (Woody 
et al. 2004), which is a recently commissioned millimeter 
array, combining the BIMA and OVRO (Owens Valley 
Radio Observatory). It consists of 6 elements of 10.4 m 
antennas and 9 elements of 6.1 m antennas. 5 In order 
to achieve a similar synthesized beam at the two wave- 
lengths, the A = 1.3 mm and 2.7 mm continuum data 
have been taken in the most compact E configuration 
and the D configuration, respectively. These two combi- 
nations of wavelengths and array configurations provide 
well matched synthesized beams, about 5" x 5". 

This moderately matched beam size at these two wave- 
lengths has not been achievable before CARMA. In in- 
tcrfcrometric observations, while high angular resolution 
can be obtained via increasing baselines of antenna ele- 
ments, there is the usual missing flux problem. This is 
because interferometric observations are only sensitive to 
size scales corresponding to the uv coverage. To mitigate 
the missing flux issue, we need either an additive single 
dish observation or well-defined uv coverage with short 
baselines. From this point of view, the most compact 
CARMA E configuration is just right to study Class 
envelope structures, since the canonical size of Class 
source envelopes is several thousands of AU correspond- 
ing to a few tens of arc-seconds in most nearby star form- 
ing regions (e.g., the Perseus molecular cloud at a dis- 
tance of 250 pc). The E configuration provides baselines 
from ~ 6 m to ~ 60 m (~ 4.6 — 46 kA at A = 1.3 mm), 
which result in a synthesized beam (angular resolution) 
of about 5" x 5". A simulation shows that our data uv 
coverage recovers fluxes well (> 50%) towards extended 
features about up to 4 times the synthesized beam size. 

CARMA has a couple of special features to realize the 
most compact E configuration. One is an anti-collision 
system installed on the 6.1m antennas, which are located 
in the inner region of the configuration. Antennas stop 
whenever they are in a danger of collision. The other 
feature is the coordinated movement. In larger configu- 
rations, D, C, and B configurations, antennas diagonally 
move (simultaneously in azimuth and elevation) to reach 
a target. However, in E configuration they go to a high 
elevation first and move in azimuth followed by a move- 

5 Recently 8 elements of 3.5 m antennas (the Sunyaev-Zel'dovich 
Array) have been merged as well. 



Grain Growth and Density Distribution in the Youngest Protostellar Systems 



3 



ment to arrive at a designated elevation, to reduce the 
collisional situations. 

The A = 2.7 mm continuum was observed in the D- 
like commissioning configuration of 2006 fall and winter 
and D configuration of 2007 summer, while the A = 1.3 
mm continuum was obtained in the E configuration of 
2007 summer. Each data set was taken with one or 
two double-side bands of a 500 MHz bandwidth in each 
single-side band for the continuum observations. Two or 
one extra bands were assigned to a CO rotational transi- 
tion (J = 2 — > 1 or J = 1 — > 0). The CO rotational tran- 
sition data are presented in another paper with other 
molecular transition data. The details of each observa- 
tion are listed in Table 1. Two and three pointing mo- 
saics have been done to better cover the larger bipolar 
outflow regions for the CO J = 2 — > 1 transition towards 
L1448 IRS 3 and L1157, respectively, at A = 1.3 mm. 
For this study, the northwest pointing data of L1448 IRS 
3 and the central pointing data of LI 157 were used. 

The Multichannel Image Reconstruction, Image Anal- 
ysis, and Display (MIRIAD, Sault et al. 1995) tools have 
been employed to reduce and analyze data. In addition to 
normal procedures (linclcngth, bandpass, flux, and gain 
calibrations), shadow-defected data have been flagged in 
the E configuration data. Shadowing indicates cases of 
an antenna's line-of-sight interrupted by other antennas 
and usually appears in low elevation observations of com- 
pact configurations. The normal effects of shadowing are 
reduction and distortion of incident antenna power and 
abnormal gain jumps. Therefore, to obtain reliable re- 
sults the shadow-defected data were flagged in the com- 
pact E configuration. 

Further special attention needs to be given on flux cal- 
ibration for studies involving flux comparison between 
different wavelengths like this study. To minimize er- 
rors caused by primary flux calibrators, we used the 
same flux calibrator (Uranus) at both wavelengths ex- 
cept L1157, which used MWC349 at A = 1.3 mm and 
Mars at A = 2.7 mm. We expect 15% and 10% uncer- 
tainties of flux calibrations at A = 1.3 mm and 2.7 mm, 
respectively, based on the CARMA commissioning task 
of flux calibration. During a commissioning period ex- 
tending to longer than 4 months, 12 calibrator (quasar) 
fluxes had been monitored by CARMA. As a result, the 
least varying case showed about 13% deviation in flux. 
When considering the intrinsic variability of quasars, it 
is expected that CARMA flux calibrations have about 
10 — 15% uncertainties. As a result, we consider 15% 
and 10% uncertainties at A = 1.3 mm and 2.7 mm, re- 
spectively. 

In addition, we make synthesized beam sizes the same 
as possible at both wavelengths, using weighting and 
tapering schemes, in order to minimize the beam size 
effect on the flux comparison. After proper weighting 
and tapering schemes, we could match the beam sizes to 
within 1%. The details of applied weighting and taper- 
ing schemes are listed in Table 2 with final synthesized 
beams. Briggs' robust parameter is used (Briggs 1995), 
which is a knob to provide intermediate weighting be- 
tween natural and uniform weighting. The parameter of 
2 gives a weighting close to natural weighting and —2 
close to uniform weighting. 



4. OBSERVATION RESULTS 

4.1. Dust opacity spectral index maps 

Total flux (F v ) of the thermal dust continuum emission 
represents the total mass (Mr) of the source, if the source 
is optically thin at the observational frequencies, 



F v « k v B v (T d ) -jp, 



(1) 



where k v , B v {Td), Mt, and D are opacity (mass absorp- 
tion coefficient) of the dust grains, blackbody radiation 
intensity of a dust temperature T d , total mass, and dis- 
tance to the source, respectively. The opacity of dust 
grains (k„) depends on dust properties such as sizes, 
components, and shapes. If the dependence is simple, 
for example a power law (k v oc v&), the dust grain prop- 
erties can be studied by observations at two frequencies. 
In addition, in the case that the Rayleigh-Jeans approx- 
imation of blackbody radiation is applicable (hv <C kT) , 
the relationship between spectral indexes of the observed 
flux densities (a) and spectral indexes of the dust grain 
opacity (0) is simply, 



F rts F 



F v ^k v B v (T d ) 



M T 

"D2" 



therefore a : 



10 + 2. 



v\P 2kT d 2 M' 



T 

7p 



(2) 



Note that this relation is valid only in the optically thin 
assumption and the Rayleigh-Jeans approximation. 

Draine (2006) showed that mainly depends on the 
size distribution of dust grains rather than their com- 
ponents and shapes; small (~ 1) is likely indicating 
dust grain size distribution up to 3A. Since our observa- 
tions are up to 3 mm, ~ 1 would suggest a grain size 
distribution up to about 1 cm. 

Figure 1 presents maps of L1448 IRS 2, L1448 IRS 3, 
and L1157. Dust continuum maps at A = 1.3 mm and 
A = 2.7 mm have been separately constructed using dif- 
ferent weightings and taperings as described in § 3 and 
Table 2 in order to have as similar synthesized beams 
as possible at the two wavelengths. Afterwards values 
of each source have been calculated using the two con- 
tinuum images. Only regions above three signal-to-noise 
ratio (SNR) levels on the both maps have been used to 
derive assuming 



_ \og{F{ Vl )/F M) 2 
log^i/^o) 



(3) 



where vi and i>q arc frequencies corresponding to 
A = 1.3 mm and A = 2.7 mm data, as listed in Table 2. 
Note that the Raylcigh- Jeans approximation and the op- 
tically thin assumption are used. In the case of an aver- 
age dust temperature of about 30 K, the upper limit of 
frequencies to which the Rayleigh-Jeans approximation 
can be applied is about 625 GHz. Since the higher fre- 
quency of our data is about 230 GHz, the assumption is 
valid for this study. However, caution should be taken 
in comparison at submillimeter wavelengths for cold 
objects such as the Class YSO envelopes. 



4 



Kwon et al. 



As shown in Figure 1, most (3 values in the three targets 
are less than 1. For a convenient comparison, the same 
gray scales have been adopted for all three maps. The 
actual ranges of (3 values are in Table 3 with the averages. 
As listed in the table, the maximum values are larger 
than 1.0. However, those large (3 values appear only on 
a few pixels of source boundaries, which may be due to 
contamination from ambient clouds. (3 and its averages 
in most regions of the three sources are similar to or less 
than 1. In the case of L 1448 IRS 3, in which three Class 
sources (3A, 3B, and 3C) exist, (3 values corresponding to 
the three sources are separately listed in Table 3. Like 
the other targets, these three sources of L1448 IRS 3 
have (3 around or less than 1. The LI 448 IRS 3 A and 3B 
fluxes are obtained simply by cutting the protuberance 
in Figure 1. Table 3 also has (3 values obtained from 
the total fluxes at the two wavelengths, which have been 
estimated in source regions limited by the three SNR 
threshold at both wavelengths. All sources except L1448 
IRS 3B have f3 values comparable to the mean values of 
the (3 maps. 

Another feature to note is that there are (3 gradients 
with radius in all sources. LI 157 has a smaller (3 in 
the northeast-to-southwest direction, roughly consistent 
with the A = 1.3 mm and 2.7 mm results of Beltran et al. 
(2004). However, it is noteworthy that they restored 
their two images with an identical beam size without 
any weighting schemes, which could cause a biased re- 
sult due to different uv coverage of the two wavelength 
data. The radial dependence of (3 is better shown in § 
4.2 and is discussed in detail for the L1448 IRS 3B case 
via modeling in § 6 

4.2. Visibility data comparison 

We have also examined (3 values in uv space, which 
is the Fourier transformed space of an image. Data 
of interferometric observations are obtained in uv space 
and called uv visibilities or just visibilities. To obtain a 
sky intensity distribution, inverse Fourier transformation 
and deconvolution (e.g., CLEANING algorithm) are em- 
ployed (e.g., Thompson et al. 2001). However, limited 
uv coverage causes difficulties, i.e., the deconvolution in- 
troduces systematic biases, especially for non-point, ex- 
tended sources. One of the best ways to overcome this 
difficulty is to investigate the visibility data in uv space 
instead. 

The results of (3 calculated in uv space are displayed 
in Figure 2. Visibilities have been vector-averaged in an- 
nuli. Since the envelope structures from our observations 
are spherical, the annulus averaging is valid. The annulus 
bin sizes are ~ 3.1 kA except when the SNR is too low, 
usually at the relatively longer baselines. This is most no- 
ticeable in L1157 at A = 2.7 mm. Although the uv cover- 
age is comparable at both wavelengths, the lower SNR at 
A = 2.7 mm requires larger bins. The (3 values are calcu- 
lated at the A = 1.3 mm bins with A = 2.7 mm visibilities 
linearly interpolated using the nearest bin values. When 
the A = 1.3 mm bin center is beyond last A = 2.7 mm bin 
center (extrapolation case), then the nearest bin value for 
A = 2.7 mm is used. 

In the case of L1448 IRS 3, only 3B is considered for 
the (3 calculation in uv space. The other two objects, 
3 A and 3C, are too small and weak to carry out the 
calculation. On the other hand, 3A and 3C should be 



removed from the visibilities to obtain the 3B data. Us- 
ing the MIRIAD task UVMODEL and image models ex- 
cluding the two components, we subtracted the 3 A and 
3C visibilities at both A = 1.3 mm and A = 2.7 mm sep- 
arately. In addition, since the A = 1.3 mm data set has 
been taken with two pointings offset from the center, we 
compensated the primary beam sensitivity loss using a 
UVMODEL multiplication. 

In Figure 2, the upper panels show amplitudes of 
A = 1.3 mm (open squares) and A = 2.7 mm cases (open 
triangles). The error bars represent the statistical stan- 
dard errors in each bin. The solid and dashed lines 
present the best fit models described in § 5 and Fig- 
ure 3. The lower panels show (3 values with uv distance, 
calculated by equation (3). The open circles indicate (3 
values calculated from the uv visibilities shown on the up- 
per panels. The error bars with caps on the open circles 
represent (3 value ranges corresponding to the statistical 
amplitude errors of the upper panels. The filled circles 
and error bars without caps show the effect that the ab- 
solute flux calibration uncertainty has on the calculation 
of (3. We adopt 15% flux calibration uncertainties for 
A = 1.3 mm data and 10% for A = 2.7 mm data, as dis- 
cussed in § 3. The larger (3 points indicate the case in 
which 15% higher fluxes at A = 1.3 mm and 10% lower 
fluxes at A = 2.7 mm arc considered and vise verse for 
the lower (3 points. The (3 ranges are around ±0.35, as 

log(1.15/0.90)/log(^i/^ ) ~ °- 35 where "i/^o ~ 2 (refer 
to eq. 3). 

Two main features should be noted in Figure 2. One 
point is that the (3 values are around 1 or less than 1 in 
all three objects. It is arguably true even when consider- 
ing the absolute flux calibration uncertainties. The other 
point is the radial dependences of (3. In L1448 IRS 2 and 
L1157, (3 arguably decreases on smaller scales (larger uv 
distances). L1448 IRS 3B, however, distinctly presents a 
radial dependence. The (3 variation is fit with the loga- 
rithmic function of (3{C,) = 1.0 — 0.57 log(C), where £ is 
the uv distance in units of kA. When assuming power-law 
distributions of density and temperature of envelopes as 
discussed in § 5, the distributions of the intensity inte- 
grated along line-of-sight as well as the radial intensity 
follow a power-law under the optically thin assumption 
and Raylcigh- Jeans approximation (Adams 1991). When 
ignoring primary beam effects of interferometers and as- 
suming infinite size envelopes, the visibilities are also in a 
power-law (e.g. Harvey et al. 2003; Looney et al. 2003). 
As (3 is obtained from equation (3) here, we assume a 
logarithmic function of (3((). There arc a few possible 
interpretations to explain this radial dependence of (3. 
It could be caused by increasing the fraction of optically 
thick emission on smaller scales due to the denser central 
region. Beckwith et al. (1990) discussed that the opti- 
cally thick emission fraction (A) decreases (3 by a factor 
of (1 + A), i.e., (3 ~ /? /(l + A). Similarly, it could be due 
to an optically thick, unresolved, deeply embedded disk 
structure at the center. On the other hand, it could indi- 
cate a faster grain growth in the denser central region or 
dust grain segregation suggested by some star formation 
theories, for example, ambipolar diffusion in magneti- 
cally supported molecular cloud (Ciolek & Mouschovias 
1996). The radial dependence is discussed in more detail 
in § 6. 
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5. MODELING IN UV SPACE 

As mentioned in § 4.2, images of extended features con- 
structed from interfcromctric observations may be biased 
due to limited uv coverage. In contrast, comparing visi- 
bility data against source models transformed to the vis- 
ibility plane (including the primary beam modification, 
Fourier transformation, and visibility sampling), is not 
prone to these imaging deconvolution biases. Therefore, 
we carry out envelope modeling in uv space rather than 
in image space. In other words, we compare observation 
visibilities with model visibilities sampled over the ob- 
servation uv coverage, after obtaining uv models by the 
Fourier-transformation of image models. 

We assume that the temperature distribution of dust 
grains is in radiative equilibrium with the central pro- 
tostar, ignoring heating by gas and cosmic rays (Spitzcr 
1978, p 193): 



poo poo 

■ \ Q a {v)u v dv = 4tt / Q a (v)B v {T d )dv, 
Jo Jo 



(4) 



where Q a (v), u v , B v (T d ), and c are absorption efficiency 
factor, radiation energy density, black body radiation in- 
tensity of temperature T d , and speed of light, respec- 
tively. The radiation energy density (u u ) at a distance 
r from the center can be expressed as ttB„(T*) Rl/r 2 , 
where T* and i?* are an effective temperature and a ra- 
dius of a central protostar. Assuming Q a (v) oc zV 3 , equa- 
tion (4) gives a temperature distribution of dust grains 
(Beckwith ct al. 1990), 



W = T.(i^) 



1 i?*\ 2/(4+/3) 



(5) 



Again, (3 is the dust grain opacity spectral index (k„ = 
ko{v/vo)P). This equation can also be formulated with a 
grain temperature T at a distance i? from the central 
protostellar luminosity L 0} as (e.g., Looney et al. 2003) 



_R Q s 2/(4+0) ,L „N 1/(4+0) 



w= T °(v) U 



(6) 



Although the inner region, which might be optically 
thick, could have a sharper temperature gradient than 
this relation (e.g., Wolfire & Cassinelli 1986; Looney et al. 
2003), it is limited at the very central region, and our 
results arc not sensitive to the possibility, as further dis- 
cussed in § 6. 

Some previous studies (e.g., Harvey et al. 2003) con- 
sidered the external heating by the interstellar radiation 
field, using a temperature lower limit of 10 K. However, 
we do not explicitly include this effect, since the temper- 
ature lower limit is uncertain and the lowest temperature 
of our modeling is comparable, about 7.3 K at r — 7000 
AU when adopting T Q = 100 K at r = 10 AU. In ad- 
dition, tests show that the temperature lower limit does 
not change our results, as previous studies also reported 
(e.g., Harvey et al. 2003). The outer envelope heated ex- 
ternally by the interstellar radiation field would be the 
main intensity component in sources without any central 
heating objects, but in Class YSOs the central high 
temperature region drives the emission. Besides, inter- 
ferometric observations are not so sensitive to the outer 
envelope, where the effect of the temperature lower limit 
is largest. 



The power-law density distribution is assumed for en- 
velopes, p(r) = po( r / r o)~ p - Therefore, the intensity of 
envelopes on the plane of the sky is calculated as 



lu = f B v {T d {r)) e- r » dr v = j B v {T d (r)) e 



p(r) k v dL, 
(7) 

where L indicates the line-of-sight from the observer and 
the optical depth r v = k v p(r) dL' . Spherical en- 
velopes with an outer radius of R out and with an inner 
hole of a radius of Ri n are assumed. Therefore, the den- 
sity distribution can be expressed with the total envelope 
mass Mt (when p ^ 3) as 

Mt — / p(r) Airr 2 dr 
J^ n 

4-7T 



3-p 
p(r)=p r p r 



(Rout 3 - p - R in 3 - p ) Por p 



47T 



(8) 
(9) 



Substituting the density expression with the total enve- 
lope mass into the optical depth of equation (7) shows 
a coupling of Mt and kq — in the case that the enve- 
lope is optically thin and the Raylcigh- Jeans approxima- 
tion {B u {T d {r)) « 2kT d (r)/\ 2 )'is valid, T is also cou- 
pled. Normally the envelopes of this stage YSO are opti- 
cally thin in the A = 1.3 mm and 2.7 mm continua except 
the very central regions (within a few tens of AU) and 
the Raylcigh- Jeans approximation is applicable, which 
means that the Mt, Ko, and To are all likely coupled. 
However, note that the optically thin assumption and 
Rayleigh- Jeans approximation, which are assumed in /3 
calculations of observational data in § 4.1 and § 4.2, are 
not assumed in the modeling to avoid biases. Here we 
just intend to point out that the three parameters are 
likely to be coupled. 

After constructing intensity image models, they are 
corrected by three different CARMA primary beams, 
which correspond to baselines of two 10.4 m antennas, 
two 6.1 m antennas, and 10.4 m and 6.1 m antennas. 
The three primary-beam corrected images are Fourier- 
transformed into uv space and model visibilities are sam- 
pled over the actual observational uv coverage of the 
three different baselines. Comparison between model and 
observation visibilities is done by vector averaged values 
in annulus bins. Although bipolar outflows at this stage 
carve a cavity (e.g., Seale & Looney 2008), the effect is 
minor (Chandler & Richer 2000), especially at our inter- 
mediate angular resolution. In addition to the bipolar 
outflow effect, envelopes might be clumpy. However, the 
effect on our modeling is also insignificant, since the an- 
gular resolution of our data is intermediate and annulus- 
averaged values are used for the comparison of models 
and data. 

Parameters involved in our modeling are p (power-law 
density index), (3 (opacity spectral index), Mt (envelope 
total mass), kq (opacity coefficient at ^o), To (grain tem- 
perature at Ro), Rin and R ou t (inner and outer radii 
of envelopes), and F pt (a central point source flux at 
A = 2.7 mm). Among these, two parameters are fixed: 
k = 0.0114 cm 2 g- 1 at v = 230 GHz and T = 100 K 
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at Rq = 10 AU. As discussed, the ko and Mt are cou- 
pled (and so To is mostly), so we cannot well constrain 
these parameters simultaneously. The T at R n corre- 
sponds to a central object luminosity of 1.67 Lq and the 
K a at v = 230 GHz is the average of = 1 and 2 cases 
in k v — 0.1 (i^/1200 GHz)^ 3 , assuming a gas-to-dust mass 
ratio of 100 (e.g., Hildebrand 1983; Beckwith et al. 1990). 
Ossenkopf & Henning (1994) also reported k « 0.01 cm 2 
g _1 at A = 1.3 mm for dense protostellar cores via dust 
coagulation model calculation, when using a gas-to-dust 
mass ratio of 100. Note that ko is not very well known 
and has a large uncertainty (e.g., Hildebrand 1983; Beck- 
with & Sargent 1991) so we need to pay attention to the 
fact that the total mass Mt could have a large uncer- 
tainty. Mt can also be scaled by the presumed T . 

The central point source flux (F pt ) is designed to sim- 
ulate an unresolved central disk structure. We assumed 
that the point sources are optically thick so that the flux 
density spectral index is 2 under the Rayleigh- Jeans ap- 
proximation, meaning = 0. In the case of L1448 IRS 
2 there is no point source required, since there is no flat 
visibility amplitude on the small scales, particularly at 
A = 1.3 mm in Figure 2. It may indicate that the central 
disk structure of the source is not so significant. In con- 
trast, L1157 has a flat profile on the small scales, which 
means a compact structure at the center. Therefore, a 
point source is adopted to fit the data. Indeed, Beltran 
et al. (2004) reported a compact component (size < 1") 
of 25 mjy and 78 mJy at A = 2.7 mm and 1.3 mm, re- 
spectively. On the other hand, the point source of L1448 
IRS 3B was applied for a different reason: to simulate a 
radial dependence of 0. As shown in Figure 2, there is 
a clear radial dependence of 0, which results in no good 
fits with a constant over all scales. It is why an opti- 
cally thick point source is considered, although there is 
no flat feature on the small scales. Note that even higher 
angular resolution observations have not detected such a 
point source signature (Looney et al. 2003). We further 
discuss the radial dependence of LI 448 IRS 3B in § 6. 

In order to find good fit models, we search grids of 
parameters, p, 0, Mt, Rin, Rout, and F pt . Parameter 
set information of the three sources is listed in Table 4. 
On each grid point of parameters, the reduced x 2 (x 2 ) 
has been calculated. The two wavelength data were used 
simultaneously for fitting. Note that the absolute xl 
values particularly in L1448 IRS 3B (~ 8.7) are large, 
compared to L1448 IRS 2 (~ 1.6) and L1157 (~ 1.5). 
This is because the relatively small standard errors due 
to the high brightness of L1448 IRS 3B make fitting very 
difficult. The L1448 IRS 3B data may have imperfect 
exclusion of the companion L1448 IRS 3A, which might 
cause a difficulty in fitting. However, it is unlikely to be 
the main effect, since the companion is relatively weak 
and we subtracted the component as mentioned in § 4.2. 
In addition, the vector averaging in annuli minimizes the 
effect. On the other hand, it may indicate that the simple 
power-law model is not appropriate to explain high SNR 
observations (e.g., Chiang et al. 2008). 

We adopt a likelihood calculation to constrain p and 0, 
instead of reporting large ranges of each parameter to fit 
the data. Reporting good fit parameter ranges could bias 
the impression of the results, since each parameter value 
in the range comes from different combinations of the 



other parameters. The likelihood function we adopt is 
cxp(— x 2 /2), since the annulus averaged visibilities have 
a Gaussian distribution based on the central limit theo- 
rem. As we want to constrain p and 0, the likelihoods 
of all grid points with common p and are summed. 
The sum now indicates the likelihood of a set of p and 
0. Finally, it is normalized by the total sum of the like- 
lihoods in each plot of Figure 3, which means that the 
plots are comparable to probability density distributions 
of p and 0. Note that we do not consider the absolute 
flux calibration uncertainties for fitting. In other words, 
we use data points marked with open symbols in Figure 
2. Note that while systematic changes of absolute fluxes 
in the same direction at both A = 1.3 mm and 2.7 mm 
affect the total mass Mt, the opposite direction changes 
mainly influence 0. We estimate that the maximum 
ranges caused by the absolute flux calibration uncertain- 
ties are ±0.35, as mentioned in § 4.2. 

We present the most likely and p in Figure 3. As 
clearly shown in the figure, of the three sources are 
most likely to be around 1 even in the modeling with- 
out the optically thin assumption and Rayleigh- Jeans ap- 
proximation. These are the first clear modeling results 
showing the of Class YSOs. The contours in Figure 
3 indicate likelihood levels of 90% down to 10% of the 
peak in steps of 10% and the triangles and circles mark 
the p and pairs of the best fit models and likelihood 
weighted averages of individual parameters, respectively. 
Note that, therefore, the combinations of the weighted 
averages are not necessarily the best fit. Since a model 
with a point source is not the best one for L1448 IRS 3B, 
its contours are drawn in dashed lines. (The best model is 
discussed in § 6 and displayed in Figure 6.) The broader 
distribution in p of LI 157 is due to the adopted point 
sources. As having a point source implies a density gradi- 
ent, it lowers the density index. The two dotted contours 
in the LI 157 plot present 90% and 80% of the peak likeli- 
hood based on all models in the whole range of the point 
source fluxes F pt (0.000 - 0.035 Jy at A = 2.7 mm) listed 
in Table 4. In contrast, the solid contours of LI 157 in 
Figure 3 show the likelihood distribution obtained from 
models in a limited range of F pt (0.015 — 0.025 Jy) around 
the likelihood weighted average (0.019 Jy at A = 2.7 mm 
and 0.078 Jy at A = 1.3 mm), which is consistent with 
the compact component flux measured by Beltran et al. 
(2004). 

While the power-law density indexes of L1448 IRS 2 
and L1157 are around 1.8 and 1.7, respectively, that of 
L1448 IRS 3B is around 2.1. The density index of L1448 
IRS 3B is consistent with the lower limit of Looney et al. 
(2003) using BIMA data and the L1157 result is consis- 
tent with that of Looney et al. (2007) using Spitzer IRAC 
absorption features. The density distribution of L1448 
IRS 2 has not been studied. It is interesting to note that 
star formation theories have suggested density indexes 
between 1.5 and 2.0; "inside-out" collapse models (Shu 
1977) suggested 1.5 for the inside free-fall region and 2.0 
for the outside isothermal envelope, where the expansion 
wave does not reach yet, and ambipolar diffusion mod- 
els (e.g., Mouschovias 1991; Tassis & Mouschovias 2005) 
suggested around 1.7 but with the very inner regions 
dependent on magnetically controlled accretion bursts. 
Although we do not attempt to constrain the star for- 
mation theories in this paper, the difference in density 
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indexes between L1448 IRS 3B and the others is note- 
worthy. The difference even increases in the better model 
of L1448 IRS 3B in § 6. 

It is important to note that the constraints on the inner 
and outer radii are not very strong. While the inner 
radius of L1448 IRS 3B is likely to be 10 AU rather than 
20 AU, there is no likely inner radius for L1448 IRS 2 and 
L1157 in the parameter search space. In addition, while 
the outer radius of L1157 is likely around 2000 - 2500 
AU, the outer radii of L1448 IRS 2 and L1448 IRS 3B 
cannot be constrained well due to lack of sensitivity of 
the data toward large scales. We can only say that the 
preferred fits for these two sources have a larger outer 
radius. The values given in Table 4 are limited by our 
parameter search space. 

6. RADIAL DEPENDENCE OF fi 

We verify the radial dependence of (3 that is shown in 
L1448 IRS 3B and attempt a modeling with j3 as a func- 
tion of radius in this section. This result is the first evi- 
dence to clearly show a radial dependence of (3 in Class 
YSOs via uv modeling. Some previous studies have 
suggested a radial dependence of (3, for example, in dust 
cores of NGC 2024 (Visser et al. 1998), the Class source 
HH211-mm (Chandler & Richer 2000), and four Class I 
sources (Hogerheijde & Sandell 2000). However, the re- 
sults are not clear and it could be due to other effects such 
as an optical thickness effect or an improper considera- 
tion of temperature effects, since their results are based 
on submillimeter wavelength observations, in which the 
[3 evaluation is more sensitive to the temperature. 

As mentioned in § 5, an optically thick point source has 
been adopted to fit L1448 IRS 3B data. First, to verify 
that the point source should be optically thick to imply a 
radial variation of (3, we tested the case of a point source 
with the same (3 to that of its envelope. As expected, a 
point source with the same (3 as the envelope requires a 
smaller (3 to fit the data (Fig. 4). The dashed contours 
in Figure 4 are 80%, 60%, and 40% of the peak value in 
the likelihood distribution of the same models in Figure 
3 and the dotted contours are 80%, 60%, and 40% of the 
likelihood peak in the new models with a point source 
having the same (3 to the envelope. A parameter space 
of p : 1.9 — 2.4 and (3 : 0.4 — 0.9 with the other param- 
eter ranges the same as the optically thick point source 
models, except Ri n which was fixed at 10 AU, has been 
searched. In addition to the smaller (3, it is noteworthy 
that there is no "good" fit. The "best fit" gives xt ~ d, 
which is much worse than the case of the optically thick 
point source case (xt ~ 8.7). This is expected as there is 
no good way to well fit the two wavelength data simul- 
taneously without a variable (3 along radius. Note that 
the differences between the two wavelength amplitudes 
are only sensitive to [3. Since we assume a constant (3 
for the point source and the envelope in the new model, 
the differences between the two wavelength amplitudes 
along radius can be caused only by the optically thick 
emission due to the density increase of the inner enve- 
lope. As the new model is worse than the optically thick 
point source model, this test also implies that the opti- 
cally thick emission, purely due to the density increase 
of the inner envelope of L1448 IRS 3B, is not significant 
enough to explain the (3 variation in the data. 

Similarly, better (probably more "realistic" , a sharper 



temperature gradient in inner regions) temperature dis- 
tributions such as of Looney et al. (2003) and Chiang 
et al. (2008) cannot fit the data either. We tested simu- 
lated temperature distributions similar to those studies 
and verified that they do not provide radially variable 
differences between the two wavelengths. The dotted 
line in Figure 5 is an example of fitting models with the 
better temperature distribution but with a constant f3 
over radius. As shown, it does not produce the variable 
amplitude differences with radius between the two wave- 
lengths. It is understandable since the inner regions are 
hotter resulting in a valid Rayleigh-Jeans approximation, 
i.e. no slope change between the two wavelengths due to 
temperature variation. 

Finally, we construct a model to simulate the variable 
(3 as a function of radius, based on grain growth. A point 
source of L1448 IRS 3B seems to be weaker than 20 mJy 
if it existed, according to Looney et al. (2003), whose 
data went to ~ 400 kA at A = 2.7 mm. Therefore, mod- 
eling with an optically thick point source is not the best 
way, although it provides a relatively "good" fit for our 
intermediate angular resolution data. For this reason, we 
do not consider a point source in the following model. 

We assume grain growth by gas accretion onto grain 
surfaces. Grains can grow by gas accretion and coagu- 
lation and can be destroyed or denuded by grain-grain 
collisions and heating mechanisms such as cosmic rays, 
central protostellar radiation, and bipolar outflow shock 
waves (e.g., Draine 1985). To address grain growth fully, 
these growth and destruction mechanisms may need to 
be taken into account together. However, we presume 
only grain growth by gas accretion without consider- 
ing any destruction mechanisms for simplicity. Coagula- 
tion might contribute significantly in the dense envelopes 
but its efficiency is uncertain (e.g., Flower et al. 2005). 
Grain growth by coagulation requires relative grain mo- 
tion, which can be introduced by various mechanisms. 
Relative velocities caused by thermal movement, am- 
bipolar diffusion, and radiation pressure lead to grain 
coagulation rather than grain shattering; the velocities 
are smaller than the critical velocities, which are the 
upper limits of velocity for grain coagulation depend- 
ing on grain properties such as size, composition, and 
shape. The critical velocities have been studied theo- 
retically (e.g., Chokshi et al. 1993) and experimentally 
(e.g., Blum 2000; Poppe et al. 2000). However, the ve- 
locity is too small to consider coagulation as an efficient 
mechanism for grain growth (Draine 1985). Alterna- 
tively, hydrodynamically or magneto-hydrodynamically 
induced turbulence (e.g., Yan et al. 2004) could bring a 
faster relative velocity of grains. However, it depends 
on the maximum velocity at the incident scale, which 
is highly uncertain, and it may also lead to grain de- 
struction due to high velocities. In addition, even when 
considering the fastest relative velocity of grains for co- 
agulation (the critical velocity), coagulation may not be 
as efficient as gas accretion (Flower et al. 2005). 

The grain growth rate by gas accretion has a re- 
lationship with density and temperature distributions, 
da/dt oc wp oc T 1 / 2 /?, where a and w indicate a grain 
size and a colliding gas velocity (Spitzer 1978, p 208). 
Note that although we assume only grain growth by gas 
accretion, grain growth rate by coagulation has a simi- 
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lar relationship with the density and relative velocity of 
grains instead of gas density and velocity. Overall, this 
formulation is arguably valid for a general description 
of grain growth, in a well-mixed gas and dust region. 
In addition to the grain growth rate, we simply assume 
that P is inversely proportional to the maximum grain 
size (Draine 2006). Therefore, after some time period, 
P(r) is inversely proportional to the product of the den- 
sity distribution and the square root of the temperature 
distribution, 



P(r) 



Pout(r/Rp) p (T d (^)/T d (r)) 1 /2 whe re r < Rp 
Pout where r > Rp. 



We fix p out — 1.7 (c.g, Draine 2006) and instead in- 
troduce R[3 for an adjustment of the radial dependence. 
In addition, we allow the temperature distribution to 
change along p(r). However, T d (r) = T (i? /r) 2/(4+/3(r)) 
is not a monotonic function, i.e., presumably not real- 
istic. Therefore, we design a temperature distribution 
smoothly changing from a case of p — to a case of 
P = Pout around Rp, 



Wi(r)T 1 {r) + W 2 {r)T 2 {r) 
W 1 (r)+W 2 (r) 



(10) 



where Wi{r) = Rp/r, W 2 (r) = r/Rp, Ti(r) = 
T (i? /r) 2 / 4 , and T 2 (r) = T (i? /r) 2 /( 4 +^"'). We rec- 
ognize that the temperature distribution might not be 
the best one corresponding to the variable p. However, 
we point out that the temperature distribution mainly 
changes the flux density profiles, not the differences be- 
tween flux densities of the two wavelengths (Fig. 5). 
Therefore, the modeling here focusing on the variable p, 
which is implied for the variable differences of the flux 
densities along radius, is not sensitive to the temper- 
ature distribution. We searched a parameter space of 
p, Mt, Rout, and Rp with the other fixed parameters 
(p out = 1.7, R in = 10 AU, F pt = 0.0 Jy, T - 100 K 
at i? = 10 AU) as listed in Table 5. Figure 6 shows 
the result, a likelihood distribution on p vs. Rp. The 
p and Rp are most likely to be about 2.6 and 400 AU, 
respectively. The parameter set of the best fit model 
(xl ~ 7.1) is p = 2.6, M T = 2.20 M Q , Rp = 400 AU, 
and Rout — 4500 AU and the averages weighted by the 
likelihood are p = 2.59, M T = 2.51 M Q , Rp = 420 AU, 
and Rout — 5900 AU. The best fit model is plotted in 
Figure 5 overlaid with the observational data. 

In this model, the best fit suggests an envelope that 
is mostly "interstellar medium grains" (small grains 
with P ~ 1.7), with grain growth at the very center, 



Rt 



< 



400 AU, which is approximately the smallest 



structure sensitivity of these observations. It is impor- 
tant to note that this is not equivalent to models of an 
"interstellar medium grain" envelope with a point source 
of a smaller P value, as those models do not fit (Fig. 4), 
and in addition, such a bright point source at A = 2.7 mm 
is not consistent with the results of Looney et al. (2003). 

The p value (~ 2.6) is larger than the value (~ 2.1) 
obtained in § 5 assuming an optically thick point source. 
This is understandable because applying a point source 
itself causes a density gradient, as mentioned in § 5 for 
L1157. Actually, this p value is more consistent with 
the results of Looney et al. (2003) using larger uv cover- 



age data and a higher angular resolution at A = 2.7 mm. 
Based on the facts that the data of L1448 IRS 3B do not 
have a point source feature and that this model has a 
smaller xl ~ 7.1, we argue that the larger p from this 
model is more reliable. 

To understand the large difference between p values of 
L1448 IRS 3B and the other two sources, we focus on the 
differences of the apparent properties. While L1448 IRS 
2 and L1157 are isolated and have a very large bipolar 
outflow (~ 5'), L1448 IRS 3B is in a "binary system" 
and its bipolar outflow is not so extended (e.g., Kwon 
et al. 2006). These facts imply that the density distri- 
bution could be steeper in binary and/or younger (based 
on the kinematic time scales of bipolar outflows) YSOs 
such as L1448 IRS 3B. Looney et al. (2003), who have 
carried out uv modeling towards 6 sources, have also re- 
ported relatively steeper density distributions for bright 
YSOs of "binary systems" such as NGC 1333 IRAS 4B 
and L1448 IRS 3B. However, density indexes larger than 
2 are somewhat puzzling, since they indicate expansion 
rather than collapse, i.e., the thermal pressure gradient 
exceeds the gravitational force. However, we might be 
able to connect this aspect to their binarity, in which 
the outer envelope is affected by the companion, or their 
youngncss, in which the envelope is affected by the bipo- 
lar outflow momentum. Detailed theoretical studies are 
needed to understand this. 

The Rp value indicates an outer limit where grain 
growth mainly occurs. According to Spitzer (1978), the 
grain growth rate by gas accretion in the diffuse inter- 
stellar medium (T = 80 K, n# = 20 cm~ 3 ) is given by, 



da 



2 x 10 



12 



T 1\V2 



80 K n 



1 20 cm" 3 J year' k ' 



assuming a typical dielectric grain density and a cos- 
mic composition gas. The £ a is a sticking probability, 
and the ^ is the mean gas particle weight. Although 
grain growth in dense regions such as the central regions 
of Class YSO envelopes could be different, it is ap- 
plicable as discussed before. Simply compensating for 
our temperature (<~ 40 K), the mean gas particle weight 
increase (two-atomic molecular gas rather than atomic 
gas), and density (n H ~ 10 9 cm" 3 at 200 AU), we can 
obtain da/dt = 5 x 10~ 5 £ a (mm/year). When accepting 
£ a = l, 6 this implies that a time scale of 10 4 years, com- 
parable to the kinematic time scales of bipolar outflows 
of Class YSOs (e.g., Bachiller et al. 2001), can result 
in about mm-size grains. Although grain growth could 
also occur in previous stages, it is much more efficient 
in the higher densities of the Class stage. Another in- 
teresting point is that less massive (i.e., less bright) and 
less steep density distribution envelopes such as those of 
LI 448 IRS 2 and LI 157 would have smaller radial regions 
for the grain growth within the same time scale. Then, 
in such sources, the variation of P may not be distinct 
nor distinguishable from a point source, as shown in § 
4.2. 

We interpreted the radial dependence of P based on 
grain growth above. However, there could be another 

6 Although Spitzer (1978) assumed £ a = 0.1 for the diffuse in- 
terstellar medium, £ a = 1 is arguably a better assumption for the 
cold and dense inner envelope regions (e.g., Flower et al. 2005). 
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effect, grain segregation. Ciolck & Mouschovias (1996) 
showed that magnetic fields in protostellar cores reduce 
abundances of small grains in the cores by a factor of its 
initial mass-to-magnetic field flux ratio. In other words, 
a stronger magnetic field with respect to the mass of 
a core causes more effective segregation. Although this 
segregation occurs while the ambipolar diffusion appears, 
before dynamical collapse, the signature footprint could 
remain in the envelopes of Class YSOs. On the other 
hand, although this effect would be minor to the features 
we have discussed because the segregation is effective to 
relatively small grains (a < 10 cm), it is noteworthy 
that it would set the initial grain distribution of Class 
YSO envelopes for more efficient growth in the central 
region. 

7. CONCLUSION 

We carried out interferometric observations towards 
three Class YSOs (L1448 IRS 2, L1448 IRS 3, and 
L1157) at A = 1.3 mm and 2.7 mm continuum using 
CARMA. The continuum at these millimeter wave- 
lengths is mainly thermal dust emission of their en- 
velopes. Our observations have been designed particu- 
larly to cover comparable uv ranges at the two wave- 
lengths, which allowed us to tackle dust grain opacity 
spectral indexes (f3) of Class YSOs, using unprece- 
dented compact configuration and high image fidelity 
Through simultaneous modeling of the two wavelength 
visibilities as well as comparisons of the images and visi- 
bilities for the first time, we found not only the j3 of Class 
YSOs but also its radial dependence. In addition, we 
addressed the single power-law density index p of Class 
YSO envelopes. 

1. We found that the dust opacity spectral index (3 of 
the earliest YSOs, so-called Class 0, is around 1. This 
implies that dust grains have significantly grown already 
at the earliest stage. 



2. Wc obtained the power-law density index p of <~ 
1.8, - 2.6, and - 1.7 for L1448 IRS 2, L1448 IRS 3B, 
and LI 157, respectively Although we did not attempt 
to constrain star formation theories, we pointed out the 
difference between that of L1448 IRS 3B and those of 
the other two. Based on different properties of L1448 
IRS 3B from the other two sources, we suggested that 
"binary system" YSOs and/or younger YSOs in terms 
of kinematic time scales of their bipolar outflows would 
have steeper density distributions. 

3. We found radial dependences of (5. In particular, 
the dependence is distinct in L1448 IRS 3B. We verified 
it by models employing (3 as a function of radius. In ad- 
dition, we discussed that the grain growth causing the 
dependence can be achieved in a time scale of 10 4 years, 
corresponding to the kinematic time scale of bipolar out- 
flows of Class YSOs. 
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TABLE 1 

Targets and Observations 
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rv ( T9finn c\\ 












Wavclc ugt li 




Plux cal. 
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r lux 
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.DCdJll blZC I\ ) 


— r T7TT» TTTQ o 
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i Qfl /ic; 1Q on 










1.3 mm 


9fl07 Ana 91 


T T I'd I'll 1 C 

\J I anno 


3C84 


4.0 


E 










0237+288 


1.2 




2.7 mm 


2006 Sep. 02 


Uranus 


0237+288 


1.6 


Comm. 11 


4" 8 X 4'.' 3 (-74°) 




2006 Sep. 12 


Uranus 


0237+288 


1.6 


Comm. 


L1448 IRS 3 


03 25 36.339 


+30 45 14.94 










1.3 mm 


2007 Aug. 19 


Uranus 


3C84 


3.9 


E 


5'.'0 x 4'.'3 (71°) 








0237+288 


1.2 




2.7 mm 


2006 Dec. 03 


Uranus 


0237+288 


1.88 


Comm. 


5'.'0 x 4'.'5 (43°) 


L1157 


20 39 06.200 


+68 02 15.90 










1.3 mm 


2007 Aug. 20 


MWC349 C 


1927+739 


0.95 


E 


4'.'6 x 3'.'8 (24°) 


2.7 mm 


2007 Jul. 12 


Mars 


1927+739 


1.6 


D 


7'.'0 x 5C6 (7°) 



a The synthesized beam in the case of natural weighting. 

b An array configuration for commissioning tasks, similar to D. Note that only part of the array was available in some cases. 
c The flux is assumed as 1.8 Jy, based on periodic CARMA flux calibrator measurements. 



TABLE 2 

Weighting and Tapering Schemes and Final Synthesized Beams 



Source 


Frequency a 
(GHz) 


Weighting 
(Robust factor) b 


Tapering Beam Size (PA) C 


Beam Ratio 
(1 mm / 3 mm) 


L1448 IRS 2 


228.60 
112.94 


0.8 
natural 


4'.'986 x 4'.'168 (-78.91°) 
4'.'826 x 4'.'277 (-74.06°) 


1.007 


L1448 IRS 3 


228.60 
112.84 


natural 
1.1 


5'.'049 x 4'.'299 (70.87°) 
4'.'951 x 4'.'412 (43.29°) 


0.994 


L1157 


228.60 
113.00 


natural 
0.0 


5'.'6 x 6'.'1 5'.'597 x 5'.'026 (-10.95°) 
5'/ 644 x 5'.'015 (-3.850°) 


0.994 



a The frequencies used for (i calculation. Refer to eq. (3). 
b Briggs' robust weighting factor (Briggs 1995). 

c Bcam size uncertainties are order of 0.1", and the values shown are to illustrate the beam size ratios. 



TABLE 3 

/3 values of the sources 





Fluxes (Jy) 






(3 maps 




Sources 


1.3 mm 


2.7 mm 




Minimum 


Maximum 


Average 


L1448 IRS 2 


0.20 


0.025 


0.95 


0.70 


1.6 


1.0 


L1448 IRS 3 












0.60 


L1448 IRS 3A 


0.090 


0.012 


0.85 


0.32 


1.7 


0.90 


L1448 IRS 3B 


1.0 


0.19 


0.35 


-0.14 a 


1.7 


0.53 


L1448 IRS 3C 


0.15 


0.026 


0.48 


0.12 


2.1 


0.59 


L1157 


0.29 


0.050 


0.49 


-0.008 a 


1.3 


0.47 



a The negative fi values are due to a bias introduced in deconvolution. 
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TABLE 4 

Model parameter sets for the three sources 



Targets 




P 


(i 


M T 




Rout 


TTi a 
■Fpt 


T b 










(M ) 


(AU) 


(AU) 


(Jy) 


(K) 


L1448 IRS 2 


A c 


1.5 - 2.0 


0.5 - 1.5 


1.00 - 2.00 
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100 




S d 


0.1 
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10 


500 


e 






best f 
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5500 
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100 




mean g 


1.79 
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1.36 


20 


5300 





100 
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0.7- 1.3 
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100 
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0.1 


0.05 


10 


500 


0.01 
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2.2 


1.1 


3.25 


10 


6500 
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100 




mean 


2.14 


0.96 


3.68 


14 


6300 


0.099 


100 


L1157 


A 


1.5-2.0 


0.5 - 1.5 


0.30 - 1.00 


10-30 


1000 - 3000 


0.000 - 0.035 


100 




8 


0.1 


0.1 


0.05 


10 


500 


0.005 






best 


1.8 


0.8 


0.55 


30 
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0.015 


100 




mean 
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0.91 


0.59 
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0.019 


100 




A h 


1.5 - 2.0 


0.5-1.5 


0.30 - 1.00 


10-30 


1000 - 3000 


0.015-0.025 


100 




mean 11 


1.72 


0.91 


0.59 


20 


2300 


0.020 


100 



a A central point source flux at A = 2.7 mm. Here the point sources are assumed as optically thick indicating f3 = 0. 
b Temperature at Rq = 10 AU 
c Parameter range searched 
d Parameter steps 
e Fixed parameter 

'Best fitting parameter set with the smallest 

g Mean of parameters weighted by the likelihood, exp(— xt/2) 

h These two lines present the cases of models with a limited point source flux range. Refer to the text. 



TABLE 5 



Model 


PARAMETER SETS 


WITH (5 AS A 


FUNCTION OF ] 
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100 




mean 
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10 


5900 


0.000 


100 
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Fig. 1. — Dust continuum in A = 1.3 mm and 2.7 mm and dust opacity spectral index (0) maps of L1448 IRS 2, L1448 IRS 3, and L1157. 
Note that most f} values are less than 1. The gray scales are the same in all three maps, although the values are distributed in different 
ranges. The statistics of the values are in Tab. 3. Synthesized beams are the same in all three maps towards each target (Tab. 2) and 
plotted on the bottom right of the (3 maps. The L1448 IRS 3A, 3B, and 3C positions came from Looney et al. (2000) and the A = 1.3 mm 
map of L1448 IRS 3 was re-centered; the pointing center was (RA, Dec) RJ (—4", 8"). The contours of dust continuum maps arc 3, 5, 9, 
17, 33, and 65 times a = 3.4 and 1.1 mjy beam -1 (A = 1.3 mm and 2.7 mm maps of L1448 IRS 2), 10 and 1.6 mjy beam -1 (A = 1.3 mm 
and 2.7 mm maps of L1448 IRS 3), and 13 and 2.4 mjy beam -1 (A = 1.3 mm and 2.7 mm maps of L1157). 
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Fig. 2. — Amplitude (upper panels) and dust opacity spectral index f3 plots (lower panels) of the three targets, L1448 IRS 2, L1448 IRS 
3B, and L1157, along uv distance. The open squares present A = 1.3 mm data and the open triangles are for A = 2.7 mm data. The error 
bars in the amplitude plots are statistical standard errors of visibilities in each bin. The solid and dashed lines present the best fit models 
described in § 5 and Fig. 3. The open circles and error bars with caps in the f3 plots indicate f3 values and distribution regions corresponding 
to the amplitude statistical errors. The filled circles and error bars without caps present cases assuming 15% higher amplitudes at A = 1.3 
mm and 10% lower amplitudes at A = 2.7 mm (resulting in the largest (3 case within absolute flux calibration uncertainty) and 15% lower 
at A = 1.3 mm and 10% higher at A = 2.7 mm (resulting in the smallest fi case). The (3 values are calculated at the uv distance bin centers 
of the A = 1.3 mm data. The visibilities of A = 2.7 mm at the positions are interpolated linearly using nearest bin values and in the case 
of extrapolation the nearest bin values are assumed. The solid line in the (3 plot of L1448 IRS 3B is a logarithmic fit to the data. Refer to 
the text for further details. 
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Fig. 3. — Model fitting results of three Class sources in likelihood. The contour levels are from 90% of the peak value with steps of 
10%. The triangles mark best fit p and f3 pairs and circles indicate likelihood weighted averages of p and (5. To indicate that the model of 
L1448 IRS 3B is not the most reliable one in this paper (refer to § 6), its contours are presented by dashed lines. The two dotted contours 
of L1157 indicate 90% and 80% of the peak likelihood based on all models in the parameter ranges of Tab. 4 and the solid contours are the 
likelihood distribution obtained from better limited models. Refer to the text for details. In the two cases, the best fit model is identical 
and the likelihood weighted averages are slightly different. 



Grain Growth and Density Distribution in the Youngest Protostellar Systems 

L1448 IRS 3B 



15 



2.4 



2.2 



■j 



I: ■/ 

•p \; v ( : 
■-. ■ \ ■■ v 



0.6 



i \ 
i \ 
O i \ 



0.8 



1 .2 



Fig. 4. — Likelihood plots of two cases, (a) dashed contours: model of a black body (optically thick) central point source same in Fig. 
3 and (b) dotted contours: model of a central point source with a j3 same as the envelope. Triangles presents best fit values and circles 
indicate likelihood weighted average values. Note that the best model of case (b) gives a worse fit (xi ~ H) than case (a) {xt ~ 8.7). The 
contour levels are 80%, 60%, and 40% of each likelihood peak value. 
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Fig. 5. — Examples of fitting models to emphasize a radial dependence of j3. The solid lines are the best fit model with (3 as a function of 
radius (xi ~ 7.1) and the dotted lines present an example of fitting models with a constant fi. Parameter sets for the best fit model (solid 
line): p = 2.6, M T = 2.20 M Q , out = 1.7, Rfi = 400 AU, R in = 10 AU, R out = 4500 AU, F pt = 0.0 Jy, T = 100 K at R = 10 AU and 
for the other one of a constant j3 (dotted line): p = 2.5, M T = 2.80 M©, = 1.0, R m = 10 AU, Rout = 4500 AU, F pt = 0.0 Jy, T = 100 
K at Ro = 10 AU. Note that the data points are the same as in Fig. 2 and the error bars are statistical standard errors. No absolute flux 
calibration uncertainties arc shown. 




Fig. 6. — Likelihood plot for models with variable fi along the envelope radius. Hp is the radius where fi = 1.7 outward. Refer to the 
text for detailed discussions. 



